Landsat SAV

Code

import numpy as np, pandas as pd, geopandas as gpd, sklearn.metrics as metrics
import geemap, ee , rasterio, sankee, warnings, pickle
from matplotlib.lines import Line2D
from rasterio import features
from geopandas import GeoDataFrame
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
from rasterio.mask import mask
from osgeo import gdal, osr, gdalconst
import os, glob, subprocess, json, fiona
import numpy.ma as ma
import rasterio, pickle
import matplotlib as mpl
import matplotlib.pyplot as plt, plotly.express as px
from sklearn import ensemble
from sklearn.cluster import DBSCAN, AgglomerativeClustering
from matplotlib_scalebar.scalebar import ScaleBar
from scipy import stats
import seaborn as sns
import statsmodels
import glob
from glob import glob
import rasterio.plot as rplot
from scipy.cluster.hierarchy import dendrogram, linkage
from shapely.geometry import Polygon, Point, box
from sklearn.metrics import confusion_matrix, cohen_kappa_score
from rasterio.plot import show
import rioxarray as rxr
from rio_cogeo.profiles import cog_profiles
from rasterio.merge import merge
from pyproj import Transformer
from sklearn.metrics import silhouette_samples, silhouette_score, ConfusionMatrixDisplay
from rasterstats import zonal_stats
from geocube.api.core import make_geocube
from glob import glob
from sklearn.cluster import KMeans, MiniBatchKMeans
warnings.filterwarnings('ignore')

Code

#ee.Authenticate()

Code

ee.Initialize()

Code

######################## subset L8 collection   ########################################
##################################################################################
def getFactorImg(factorNames, image):
    factorList = image.toDictionary().select(factorNames).values()
    return ee.Image.constant(factorList)

######################## pandas normalise   ########################################
##################################################################################
def NormalizeData(data):
    return (data - np.min(data)) / (np.max(data) - np.min(data))

######################## L8 cloud mask   ########################################
##################################################################################
def getQABits(image, start, end, mask):
    pattern = 0
    for i in range(start,end-1):
        pattern += 2**i
    return image.select([0], [mask]).bitwiseAnd(pattern).rightShift(start)


def maskQuality(image):
    QA = image.select('QA_PIXEL')
    shadowMask = getQABits(QA,3,3,'cloud_shadow')
    cloudMask = getQABits(QA,5,5,'cloud')
    cirrusMask = getQABits(QA,9,9,'cirrus_detected')
    scaleImg = getFactorImg(['REFLECTANCE_MULT_BAND_.|TEMPERATURE_MULT_BAND_ST_B10'], image)
    offsetImg = getFactorImg(['REFLECTANCE_ADD_BAND_.|TEMPERATURE_ADD_BAND_ST_B10'], image)
    scaled = image.select('SR_B.|ST_B10').multiply(scaleImg).add(offsetImg)
    
    return (image.addBands(scaled, None, True)
            .updateMask(cirrusMask)
            # .updateMask(cloudMask)
            # .updateMask(shadowMask)
            )



def prepSrL8(image):
    dilateMask = image.select('QA_PIXEL').bitwiseAnd(int('11111', 0)).eq(0)
    cirrusMask = image.select('QA_PIXEL').bitwiseAnd(int('11111', 2)).eq(0)
    cloudMask = image.select('QA_PIXEL').bitwiseAnd(int('11111', 3)).eq(0)
    shadowMask =  image.select('QA_PIXEL').bitwiseAnd(int('11111', 4)).eq(0)
    snowMask = image.select('QA_PIXEL').bitwiseAnd(int('11111', 5)).eq(0)
    saturationMask = image.select('QA_RADSAT').eq(0)
    scaleImg = getFactorImg(['REFLECTANCE_MULT_BAND_.|TEMPERATURE_MULT_BAND_ST_B10'], image)
    offsetImg = getFactorImg(['REFLECTANCE_ADD_BAND_.|TEMPERATURE_ADD_BAND_ST_B10'], image)
    scaled = image.select('SR_B.|ST_B10').multiply(scaleImg).add(offsetImg)

    return (image.addBands(scaled, None, True)
                 .updateMask(cirrusMask)
                 # .updateMask(dilateMask)
                 # .updateMask(cloudMask)
                 # .updateMask(shadowMask)
                 # .updateMask(snowMask)
                 #.updateMask(saturationMask)
           )

######################## L8 scale factors   ####################################
##################################################################################
def apply_scale_factors(image):
    opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2)
    thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0)
    return image.addBands(opticalBands, None, True).addBands(thermalBands, None, True)

######################## convert L8 to SST   ####################################
##################################################################################
def sst(img):
    thermalBands = img.select('ST_TRAD')
    ln = k2_im.divide((k1_im.divide(thermalBands).add(1).log())).multiply(0.001).clip(site).rename("SST")
    #.convolve(kern)
    return img.addBands(ln)

######################## scale image   ####################################
##################################################################################
def scale(image):
    return ee.Image(image).multiply(0.001)

def addNDAVI(image):
    ndvi = image.normalizedDifference(['SR_B1', 'SR_B5']).rename('NDAVI')
    return image.addBands(ndvi)

def addNDWI(image):
    ndvi = image.normalizedDifference(['SR_B5', 'SR_B6']).rename('NDWI')
    return image.addBands(ndvi)


# Function to add raster to map
def addras(Map, layer, mi, ma, name):
    Map.addLayer(layer,
                 {'bands': ['SR_B4', 'SR_B3', 'SR_B2'],'max': ma, 'min' : mi},
                 name=name)
    
# Function to add raster to map
def addl5(Map, layer, mi, ma, name):
    Map.addLayer(layer,
                 {'bands': ['SR_B3', 'SR_B2', 'SR_B1'],'max': ma, 'min' : mi},
                 name=name)
    
def fmask(image):
    # see https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC09_C02_T1_L2
    # Bit 0 - Fill
    # Bit 1 - Dilated Cloud
    # Bit 2 - Cirrus
    # Bit 3 - Cloud
    # Bit 4 - Cloud Shadow
    qaMask = image.select("QA_PIXEL").bitwiseAnd(int("11111", 2)).eq(0)

    # Apply the scaling factors to the appropriate bands.
    opticalBands = image.select("SR_B.").multiply(0.0000275).add(-0.2)

    # Replace the original bands with the scaled ones and apply the masks.
    return image.addBands(opticalBands, None, True).updateMask(qaMask)

# def addkNDAVI(image):
#     #Compute D2 a rename it to d2
#     var D2 = nir.subtract(red).pow(2)
#     .select([0],['d2']);
#     // Gamma, defined as 1/sigmaˆ2
#     var gamma = ee.Number(4e6).multiply(-2.0);
#     // Compute kernel (k) and KNDVI
#     var k = D2.divide(gamma).exp();
#     var kndvi = ee.Image.constant(1)
#     .subtract(k).divide(
#     ee.Image.constant(1).add(k))
#     .select([0],['knd']);
#     return image.addBands(kndvi);

Code

hab = geemap.shp_to_ee(os.path.join(os.getcwd(),
          'Data\\BOSS\\Cleaned\\BOSS_clip_cluster.shp'))

sites = geemap.shp_to_ee(os.path.join(os.getcwd(),
          'Data\\ICoAST sites\\WGS84\\ICOAST_midwest.shp'))

kal = geemap.shp_to_ee(os.path.join(os.getcwd(),
          'Data\\ICoAST sites\\Sites\\Kalbarri.shp'))

cliff = geemap.shp_to_ee(os.path.join(os.getcwd(),
          'Data\\ICoAST sites\\Sites\\Cliff_head.shp'))

mask = geemap.shp_to_ee(os.path.join(os.getcwd(),
          'Data\\Lidar mask\\Lidar_mask.shp'))

lidar_url = 'gs://cog-bucket/lidar_cog.tif'
lidar = geemap.load_GeoTIFF(lidar_url).rename('lidar')
tr = geemap.shp_to_ee(os.path.join(os.getcwd(),
          'Data\\ICoAST sites\\Sites\\Two_rocks.shp'))
jur =  geemap.shp_to_ee(os.path.join(os.getcwd(),
          'Data\\ICoAST sites\\Sites\\Jurien.shp'))

Code

start_date = '2022-01-01'
end_date = '2022-12-30'

L8_2022 = (
    ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2021-01-01'
end_date = '2021-12-30'

L8_2021 = (
    ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2020-01-01'
end_date = '2020-12-30'

L8_2020 = (
    ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2019-01-01'
end_date = '2019-12-30'

L8_2019 = (
    ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2018-01-01'
end_date = '2018-12-30'

L8_2018 = (
    ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)


start_date = '2017-01-01'
end_date = '2017-12-30'

L8_2017 = (
    ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)


start_date = '2016-01-01'
end_date = '2016-12-30'

L8_2016 = (
    ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2015-01-01'
end_date = '2015-12-30'

L8_2015 = (
    ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2014-01-01'
end_date = '2014-12-30'

L8_2014 = (
    ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2013-01-01'
end_date = '2013-12-30'

L8_2013 = (
    ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)


start_date = '2012-01-01'
end_date = '2012-12-30'

L5_2012 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2011-01-01'
end_date = '2011-12-30'

L5_2011 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2010-01-01'
end_date = '2010-12-30'

L5_2010 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2009-01-01'
end_date = '2009-12-30'

L5_2009 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2008-01-01'
end_date = '2008-12-30'

L5_2008 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2007-01-01'
end_date = '2007-12-30'

L5_2007 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2006-01-01'
end_date = '2006-12-30'

L5_2006 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2005-01-01'
end_date = '2005-12-30'

L5_2005 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

start_date = '2004-01-01'
end_date = '2004-12-30'

L5_2004 = (
    ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .filterBounds(sites)
    .filterMetadata('CLOUD_COVER', 'less_than', 5)
    .filterDate(start_date, end_date)
    .map(fmask)
    .median()
    .clip(sites)
)

Code

Map = geemap.Map(toolbar_ctrl=True, layer_ctrl=True, add_google_map = False)
#legend = LegendControl({}, name="Landsat 8 2013-2022", position="topleft")
#Map.add_control(legend)
addras(Map, L8_2022, 0, 0.05, '2022')
addras(Map, L8_2021, 0, 0.05, '2021')
addras(Map, L8_2020, 0, 0.05, '2020')
addras(Map, L8_2019, 0, 0.05, '2019')
addras(Map, L8_2018, 0, 0.05, '2018')
addras(Map, L8_2017, 0, 0.05, '2017')
addras(Map, L8_2016, 0, 0.05, '2016')
addras(Map, L8_2015, 0, 0.05, '2015')
addras(Map, L8_2014, 0, 0.05, '2014')
addras(Map, L8_2013, 0, 0.05, '2013')
addl5(Map, L5_2012, 0, 0.05, '2012')
addl5(Map, L5_2011, 0, 0.05, '2011')
addl5(Map, L5_2010, 0, 0.05, '2010')
addl5(Map, L5_2009, 0, 0.05, '2009')
addl5(Map, L5_2008, 0, 0.05, '2008')
addl5(Map, L5_2007, 0, 0.05, '2007')
addl5(Map, L5_2006, 0, 0.05, '2006')
addl5(Map, L5_2005, 0, 0.05, '2005')
addl5(Map, L5_2004, 0, 0.05, '2004')

Code

Map.centerObject(sites.geometry(), 5)
Map